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Abstract Pituitary cells of the anterior pituitary gland secrete hormones in response 
to patterns of electrical activity. Several types of pituitary cells produce short bursts 
of electrical activity which are more effective than single spikes in evoking hormone 
release. These bursts, called pseudo-plateau bursts, are unlike bursts studied mathe- 
matically in neurons (plateau bursting) and the standard fast-slow analysis used for 
plateau bursting is of limited use. Using an alternative fast-slow analysis, with one 
fast and two slow variables, we show that pseudo-plateau bursting is a canard-induced 
mixed mode oscillation. Using this technique, it is possible to determine the region of 
parameter space where bursting occurs as well as salient properties of the burst such 
as the number of spikes in the burst. The information gained from this one-fast/two- 
slow decomposition complements the information obtained from a two-fast/one- slow 
decomposition. 



WTeka 

Department of Mathematics, Florida State University, Tallahassee, FL 32306, USA 
e-mail: wteka@math.fsu.edu 

J Tabak 

Department of Biological Science, Florida State University, Tallahassee, FL 32306, USA 
e-mail: joel@neuro.fsu.edu 

T Vo • M Wechselberger 

School of Mathematics and Statistics, University of Sydney, Sydney, NSW 2006, Australia 
T Vo 

e-mail: thvo4579@mail.usyd.edu.au 

M Wechselberger 

e-mail: wm@maths.usyd.edu.au 

R Bertram (IS!) 

Department of Mathematics, and Programs in Neuroscience and Molecular Biophysics, Florida State 
University, Tallahassee, FL 32306, USA 
e-mail: bertram@math.fsu.edu 



Springer 



Page 2 of 23 



Teka et al. 



Keywords Bursting • mixed mode oscillations • folded node singularity • canards • 
mathematical model 



1 Introduction 

Bursting is a common pattern of electrical activity in excitable cells such as neurons 
and many endocrine cells. Bursting oscillations are characterized by the alternation 
between periods of fast spiking (the active phase) and quiescent periods (the silent 
phase), and accompanied by slow variations in one or more slowly changing vari- 
ables, such as the intracellular calcium concentration. Bursts are often more efficient 
than periodic spiking in evoking the release of neurotransmitter or hormone [1-3]. 

The endocrine cells of the anterior pituitary gland display bursting patterns with 
small spikes arising from a depolarized voltage [2-5]. Similar patterns have been 
observed in single pancreatic ft -cells isolated from islets [6-8]. Figure 1(a) shows 
a representative example from a GH4 pituitary cell. Several mathematical models 
have been developed for this bursting type [5, 8-10]. Prior analysis showed that the 
dynamic mechanism for this type of bursting, called pseudo-plateau bursting, is sig- 
nificantly different from that of square-wave bursting (also called plateau bursting) 
which is common in neurons [11-13]. Yet this analysis did not determine the possi- 
ble number of spikes that occur during the active phase of the burst. The goal of this 
paper is to understand the dynamics underlying pseudo-plateau bursting, with a focus 
on the origin of the spikes that occur during the active phase of the oscillation. 

Minimal models for pseudo-plateau bursting can be written as 

€iV = f(V,n,c) (1.1) 
n = g(V,n) (1.2) 
c = €2 h(V,c) (1.3) 

where V is the membrane potential, n is the fraction of activated delayed rectifier K + 
channels, and c is the cytosolic free Ca 2+ concentration. The velocity functions are 
nonlinear, and €\ and €2 are parameters that may be small. 

The variables V, n and c vary on different time scales (for details, see Section 2). 
By taking advantage of time-scale separation, the system can be divided into fast and 
slow subsystems. In the standard fast/slow analysis one considers €2 ~ 0, so that V 
and n form the fast subsystem and c represents the slow subsystem. One then studies 
the dynamics of the fast subsystem with the slow variable treated as a slowly vary- 
ing parameter [12, 15-18]. This approach has been very successful for understand- 
ing plateau bursting, such as occurs in pancreatic islets [19], pre-Botzinger neurons 
of the brain stem [20], trigeminal motoneurons [21] or neonatal CA3 hippocampal 
principal neurons [14], Fig. 1(b). It has also been useful in understanding aspects of 
pseudo-plateau bursting such as resetting properties [11], how fast subsystem man- 
ifolds affect burst termination [17], and how parameter changes convert the system 
from plateau to pseudo-plateau bursting [12]. 

An alternate approach, which we use here, is to consider €\ ~ 0, so that V is 
the sole fast variable and n and c form the slow subsystem. With this approach, we 



4^ Springer 



Journal of Mathematical Neuroscience (2011) 1:12 



Page 3 of 23 



-10 r 




(a) (b) 

Fig. 1 (a) Pseudo-plateau bursting in a GH4 pituitary cell line, (b) Plateau bursting in a neonatal CA3 
hippocampal principal neuron. Reprinted with permission from [14]. 



show that the active phase of spiking arises naturally through a canard mechanism, 
due to the existence of a folded node singularity [22-25]. Also, the transition from 
continuous spiking to bursting is easily explained, as is the change in the number of 
spikes per burst with variation of conductance parameters. Thus, the one-fast/two- 
slow variable analysis provides information that is not available from the standard 
two-fast/one-slow variable analysis in the case of pseudo-plateau bursting. 



2 The mathematical model 

We use a model of the pituitary lactotroph, which produces pseudo-plateau bursting 
over a range of parameter values [10]. To achieve a minimal form, we use the model 
without A-type K + current (I a)- It includes three variables: V (membrane potential), 
n (fraction of activated delayed rectifier K + channels), and c (cytosolic free Ca 2+ 
concentration). The equations are: 

dV 

Cm ~dt = ~ (ICa + lK + lK(Ca) + lBK) (2 ' !) 
dn _ (noo(V) -n) ^ 
dt z n 

^- = -f c (aI C a+k c c) (2.3) 
dt 

where Ic a is an inward Ca 2+ current, Ik is an outward delayed rectifying K + cur- 
rent, Ik{Co) is a small-conductance Ca 2+ -activated K + current, and Ibk is a fast- 
activating large-conductance BK-type K + current. The currents in the equations 
above are: 

Ica = gCam 00 (V)(V-V Ca ) (2.4) 
lK=g K n(V-V K ) (2.5) 

lK{Ca) = gK(Ca)Soo(c)(V - V K ) (2.6) 
I B K=gBKboc(V)(V-V K ). (2.7) 

Springer 



Page 4 of 23 



Teka et al. 



Table 1 Parameter values for the lactotroph model. 



Parameter 


Value 


Description 


Cm 


5 pF 


Membrane capacitance of the cell 




2nS 


Maximum conductance of Ca^ channels 


v Ca 


50 mV 


Reversal potential for Ca^ 


v m 


-20 mV 


Voltage value at midpoint of moo 




12 mV 


Slope parameter of irioc 


§K 


4nS 


Maximum conductance of K"^ channels 


v K 


-75 mV 


Reversal potential for 


Vfi 


-5 mV 


Voltage value at midpoint of n oo 


s n 


10 mV 


Slope parameter of n^o 


x n 


43 ms 


Time constant of n 


8K(Ca) 


1.7 nS 


Maximum conductance of K(Ca) channels 


Ka 


0.5 /xM 


c at midpoint of Sqo 


8BK 


0.4 nS 


Maximum conductance of BK-type K + channels 


v b 


-20 mV 


Voltage value at midpoint of /oo 


Sb 


5.6 mV 


Slope parameter of /oo 


fc 


0.01 


Fraction of free Ca 2+ ions in cytoplasm 


a 


0.0015 /iMfC -1 


Conversion from charge to concentration 


k c 


0.16 ms -1 


Rate of Ca 2+ extrusion 



The steady state activation functions are given by: 

-l 

(2.8) 

l 

(2.9) 
(2.10) 

l 

(2.11) 

Default parameter values are given in Table 1 . 

The variables V, n and c vary on different time scales. The time constant of V 
is given by t v = C m /g To taU where g To tal = gRn + gBKb<x>(V) + gCa^ooiV) + 
gK(Ca)Soo(c). During a bursting oscillation, the minimum of gTotal is 0.483 pS and 

c c 

the maximum is 3 pS. Hence, _ o m — < xy < ■ m — , or 1.7 ms < Ty < 10.4 ms, 

^ ' max gTotal mm gTotal — v — 

for C m = 5 pF, a typical capacitance value for lactotrophs. The time constant for n is 
x n =43 ms. For the variable c, the time constant is j^- = ( Q q^q 16 ) ms = 625 ms. 
Thus, n and c change more slowly than V. This time scale separation between V and 
(c, n) can be accentuated when C m is made smaller than the default 5 pF, i.e., when 
C m -> 0, xy gets smaller and V varies much faster. Thus, we can view the capacitance 



moo(V) = y +exp 

«oo(V) = ^1 +exp 
c 2 

*oc(c) = — -y 

boc(V) = ( 1+exp 



Sin 

v n -y 



v b -V 

Sb 
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C m as a representative of the dimensionless singular perturbation parameter e\ in this 
model (Eq. 1.1). 

All numerical simulations and bifurcation diagrams (both one- and two-parameter) 
are constructed using the XPPAUT software package [26], using the Runge-Kutta in- 
tegration method, and computer codes can be downloaded from the following web- 
site: http://www.math.fsu.edu/~bertram/software/pituitary. The surface in Fig. 9 was 
constructed using the AUTO software package [27]. All graphics were produced with 
the software package MATLAB. 



3 Geometric singular perturbation theory 



3 . 1 The reduced system 

We consider the full system (Eqs. (2.1)-(2.3)) as having one fast variable V and two 
slower variables n and c. The time-scale separation can be accentuated by decreasing 
the singular perturbation parameter C m . This facilitates analysis of the system dy- 
namics [28]. In the limit C m —> 0, the trajectories of the system lie on a 2-D surface 
called the critical manifold. If we define the right hand side of Eq. (2.1) by 



f(V, c, n) = -(I Ca + Ik + I K (Ca) + Ibk) 
then the critical manifold is the surface S satisfying 

S = {(V, c, n) e R 3 : f(V, c, n) = 0}. 
The equation f(V,c,n) = 0 can be solved in explicit form for n as 



1 

n = n(c, V) = 

gK 



(^777 7T~V + 8K(Ca)Soo(c) + g#^oo( V) 

(V - Vk) 



(3.1) 



(3.2) 



. (3.3) 



The critical manifold (3.3) is a folded surface (Fig. 2) that consists of three sheets 
separated by two fold curves (L~ and L + ). The lower and upper sheets are attracting 
(fy < 0) and the middle sheet is repelling (|£ > 0). The lower (L~) and upper (L + ) 



fold curves are given by 

L ± = | (V, c, n) e R 3 : f(V, c, n) = 



0and^-(V,c,rc) = 0 

dv 



(3.4) 



This yields two constant V values and two equations for n in the form of n = n(c). 
Thus, the fold curves (L ± ) are (V ± ,c,^ ± (c)) where V~ and V + are constant V 
values. The curve L + is projected vertically (along the fast variable V) onto the lower 
sheet to obtain the projection curve P(L + ), and similarly for the (L~) projection 
onto the upper sheet. Figure 2 shows the critical manifold, the fold curves and the 
projections of the fold curves. 

The reduced flow (when C m -> 0) is described by (3.3), the differential equation 
for c (Eq. (2.3)), and a differential equation for V which can be obtained by differen- 
tiating f(V, c, n) = 0 with respect to time. That is, 



df dV _ df dc df dn 
dV~dT ~ ~dc~dt + 3n~dt 



(3.5) 
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Fig. 2 The critical manifold 
and fold curves with their 
projections for g% = 4 nS and 
g BK = 0.4 nS. The curves L~ 
and L+ are the lower and upper 
fold curves, respectively. P(L _ ) 
and P(L + ) are the projections of 
L~ and L + onto the upper and 
lower sheets of the critical 
manifold, respectively. FN is a 
folded node singularity, and SC 
(green curve) is the strong 
canard. The singular periodic 
orbit (black curve) is 
superimposed on the critical 
manifold. 




where n satisfies Eq. (3.3), and n, c satisfy Eqs. (2.2), (2.3). The two differential 
equations for the reduced system are thus 

df dv , ,df /(rtoo00-n)\a/ 

— = {-fcictlca + k c c)) — + — (3.6) 

dV dt V Jcy n dc V r„ J dn v } 

^ = -f c (aIca+k c c). (3.7) 
at 

Since |^ = 0 on L ± , the reduced system is singular along the fold curves. The sys- 
tem can be desingularized by rescaling time with r = — (^y)~ 1 t. The desingularized 
system is then 



~T = {-fcioilca + k c c)) — + — 

dx dc \ r n I an 



(3.8) 



dc df 

— = f c (aI C a+k c c)-t-. (3.9) 
dr dV 



Defining 



F(V, c, n) = (-fdalca + k c c)) ^ + ( (n °° (V) (3.10) 
v 7 ac \ T n ) an 

we have the desingularized system 

dV 

— =F(V,c,n) (3.11) 
dx 

dc df 

— = f c (alca+k c c)-!-. (3.12) 
dx dV 

The desingularized system describes the flow on the critical manifold. Because of 
the time rescaling, the flow on the middle sheet, where |y > 0, must be reversed to 
obtain the equivalent reduced flow. 
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3.2 Folded singularities and canards 

Equilibria of the desingularized system are classified as ordinary singularities and 
folded singularities. An ordinary singularity is an equilibrium of Eqs. (2.1)-(2.3) and 
satisfies 

f(V,c,n) = 0 (3.13) 
n = n QO (V) (3.14) 

c= — . (3.15) 

k c 

A folded singularity lies on a fold curve (L + or L~), and satisfies: 

f(V,c,n) = 0 (3.16) 

F(V,c,n) = 0 (3.17) 
V 

— =0. (3.18) 

dv 

A folded singularity is classified as a folded node if the eigenvalues are real and have 
the same sign, a folded saddle if the eigenvalues are real and have opposite signs, or 
a folded focus if the eigenvalues are complex [22, 23, 25, 29]. For parameter values 
used in Fig. 2, the system has a folded node (with negative eigenvalues) on L + (FN, 
blue point, in Fig. 2), and a folded focus on L~ (not shown). 

There are an infinite number of singular trajectories on the top sheet that pass 
through the folded node (FN). These are called singular canards [22]. The singular 
canard that enters the FN in the direction of the strong eigenvector is called the strong 
canard (SC, green curve, in Fig. 2). This curve and the fold curve L + delimit the 
singular funnel that consists of all initial conditions whose trajectories for the reduced 
system pass through the folded node. The singular funnel and key curves are projected 
onto the (c, V) -plane in Fig. 3. The different panels are obtained with different values 
of the parameter gK • 



3.3 Singular periodic orbits, relaxation oscillations, and mixed mode oscillations 

A singular periodic orbit (Fig. 2, black curve with arrows) can be constructed by 
solving the desingularized system for the flow on the top and bottom sheets of the 
critical manifold, and then projecting the trajectory from one sheet to the other along 
fast fibers when the trajectory reaches a fold curve. The singular periodic orbit is 
the closed curve constructed in this way. This process was discussed in detail in [22, 
28, 30]. Briefly, the trajectory moves along the bottom sheet until L~ is reached. At 
this point the reduced flow is singular (|^ =0). The quasi-steady state assumption 
/ (V,c,n) = 0 is no longer valid and there is a rapid motion away from the fold curve 
L~ . This rapid motion is seen as vertical movement to the top sheet (the dynamics 
are governed by the layer problem, see [22, 28]). The trajectory moves to a point on 
P(L~) and from there is once again governed by the desingularized equations, mov- 
ing along the top sheet until L + is reached. The fast vertical downward motion along 
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0.2 0.205 0.21 0.215 0.22 0.225 0.23 0.26 0.27 0.28 0.29 0.3 0.31 

c (pM) c (fiM) 

(a) (b) 



Fig. 3 The critical manifold is projected onto the (c, V)-plane for (a) g% = 5.1 nS, gBK — 0-4 nS and 
(b) g K =4 nS, gBK — 0-4 nS. L~ and L + are the lower and upper fold curves, respectively. P(L _ ) and 
P(L + ) are the projections of L - and L + . The shaded regions are singular funnels which are delimited by 
the curves L + and the strong canards (SC, green curves). The singular periodic orbits (black curves with 
arrows) are superimposed. FN is a folded node singularity. 8 < 0 in panel (a) and 8 > 0 in panel (b). 



fast fibers returns the trajectory to a point on P(L + ) on the bottom sheet, completing 
the cycle. 

When the singular periodic orbit reaches L~ it jumps up to a point on P(L~). If 
this point on P(L _ ) is in the singular funnel, then the orbit will move through the 
FN. Otherwise it will not. Let 8 denote the distance measured along P(L~) from the 
phase point on P(L~) of the singular periodic orbit to the strong canard (SC in Fig. 3). 
When the phase point is on the strong canard, 8 = 0. Let 8 > 0 when the phase point 
is in the singular funnel and 8 < 0 when the phase point is outside the singular funnel. 
Singular canards are produced when 8 > 0. 

In Fig. 3(a) the singular periodic orbit jumps to a point on P(L~) outside of the 
singular funnel (8 < 0), so it does not enter the FN. This orbit is a relaxation oscilla- 
tion [31]. In Fig. 3(b) 8 > 0, so the orbit is a singular canard. Away from the singular 
limit, this singular canard perturbs to an actual canard that is characterized by small 
oscillations about L + [22] . The combination of these small oscillations with the large 
oscillations that occur due to jumps between upper and lower sheets yields mixed 
mode oscillations [24, 32]. The small oscillations have zero amplitude in the singular 
case, which grows as v^m for C m sufficiently small [23]. A discriminating condition 
between relaxation and mixed mode oscillations is 8 = 0, where the singular periodic 
orbit jumps to P(L~) on the SC curve. 

When C m > 0 the full system (Eqs. (2.1)-(2.3)) produces spiking for 8 < 0 and 
mixed mode oscillations for 8 > 0. Figure 4 shows these two different cases for 
g BK = 0.4 nS. For gK = 5.1 nS (8 < 0 in Fig. 3(a)), the nearly-singular periodic 
orbit produced when C m = 0.001 pF (Fig. 4(a)) perturbs to continuous spiking when 
C m = 10 pF (Fig. 4(e)). When gx = 4 nS the singular periodic orbit enters the singu- 
lar funnel (Fig. 3(b)), so when C m is increased the singular orbit transforms to mixed 
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C = 0.001 pF 

m 1 r- 



g K = 5.1nS 



g K = 4nS 




C m = 0.5pF 5 
-15 
-35 
-55 
-75 



> 

a 



■ 



i ^ i 



"i r 



i H-^ i 



0 , v 0.5 
(e) 



Time (sec) 



2 0 0.5 1 1.5 

w Time (sec) 



0 /x 0.1 0.2 0.3 0.4 0.5 0.6 0 _ 0.1 0.2 0.3 0.4 0.5 0.6 




Fig. 4 Nearly- singular periodic orbits perturb to continuous spiking or mixed mode oscillations. In both 
cases gBK — 0-4 nS, and g% = 5.1 nS in the left column, g% = 4 nS in the right column. C m is increased 
from top row to bottom row. (a), (c), (e) The singular periodic orbit does not enter the singular funnel 
(8 < 0) so it perturbs to continuous spiking, (b), (d), (f) The singular periodic orbit enters the singular 
funnel (8 > 0) so it perturbs to mixed mode oscillations or pseudo plateau bursting. 



mode oscillations. For C m = 0.5 pF mixed mode oscillations with small spikes are 
produced (Fig. 4(d)). As C m is increased to 10 pF, mixed mode oscillations with 
larger spikes are produced. This is the genesis of pseudo-plateau bursting (Fig. 4(f)). 



4 Analysis of the desingularized system and folded nodes 

We next discuss the singularities of the desingularized system for a range of gK 
and gBK values (Fig. 5). The system (with gBK = 0.4 nS) has a single-branched V- 
nullcline (green curve) that satisfies F(V, c, n) = 0 and a three-branched c-nullcline 
(orange curves) L~, L + and CN1. The curves L~, L + satisfy |£ = 0, and are the 
same as the fold curves in Fig. 3. The curve CN1 satisfies alc a +k c = 0. There are 
folded singularities that are located at intersections of the V-nullcline with L~ or L + , 
and ordinary singularities that are located at intersections with CN1. For fixed gBK, 
changing gx affects the position of the V-nullcline but not the c-nullcline. 

For values gK < 0.5131 nS, there is a stable node on CN1 (Ai), which would be 
on the top sheet of the critical manifold. There are also two folded saddles on L + (B i 
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-0.75 -0.6 -0.45 -0.3 -0.15 0 0.15 0.3 0.45 0.6 0.75 



Fig. 5 V-nullclines (green), the three-branched c-nullcline (orange), and singularities for g^K — 0-4 nS 
and different values of g% (units in nS). Filled circles represent stable singularities and unfilled circles 
represent unstable singularities. Red circles (filled or unfilled) are ordinary singularities. Filled and unfilled 
circles in blue are folded nodes and folded saddles, respectively. Filled circles in cyan are folded foci. The 
points TR1 and TR2 are transcritical bifurcations (type II folded saddle-node bifurcations) and SN1 and 
SN2 are standard saddle-node bifurcations (type I folded saddle-node bifurcations). 

and Ci) and two folded foci on L~ (Di and Ei). When gK is increased to 0.5131 nS 
the stable node Ai moves down and to the left and the folded saddle Bi moves to the 
left. These two equilibria coalesce at a transcritical bifurcation (TR1). This transcrit- 
ical bifurcation corresponds to a bifurcation of folded singularities called a type II 
folded saddle-node [22, 30, 33]. Following this bifurcation, the folded singularity is 
a folded node. For gK = 4 nS, the equilibria on L + are the folded node (B3) and the 
folded saddle (C3). The equilibrium on CN1 (A3) is now a saddle point. There is no 
qualitative change of equilibria on L~ . 

When gK is increased to 7.588 nS the equilibria B3 and C3 coalesce at a saddle- 
node bifurcation point (SN1). This is a standard saddle-node bifurcation of folded 
singularities and is called a type I folded saddle-node [22, 30, 33]. As gK is increased 
to 43.1 nS, the folded focus D5 moves to the left and changes to a folded node at D6. 
The saddle points on CN1 move downward and to the left as gK is increased. For 
gK = 129.2 nS, the saddle point A6 coalesces with the fold node D6 at a second 
transcritical bifurcation (TR2); again a type II folded saddle-node. Beyond this, the 
ordinary singularity (As, A9) is stable and the folded singularity becomes a folded 
saddle. Moreover the folded focus E6 has become a folded node (E7). As gK is in- 
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creased further to 137.2 nS, there is a second type I saddle-node bifurcation (SN2) at 
which the folded node and the folded saddle coalesce and disappear. For the values 
gK > 137.2 nS, the only equilibrium is on CN1 and is an ordinary stable node (A9). 
This is on the bottom sheet of the critical manifold. 

Varying gBK slightly affects the V-nullcline and strongly affects the c-nullcline 
in the (c, V)-phase plane. Increasing gBK moves the fold curves together, eventually 
taking the fold out of the critical manifold. Figure 6 shows qualitative changes in the 
equilibria when gBK is varied, with gx = 7.588 nS. When gBK =0.2 nS there is 
a saddle point on CN1 (A) and two folded foci (D and E) on L~ (Fig. 6(a)). When 
gBK is increased to 0.4 nS, the curve L + moves down and a type I folded saddle-node 
bifurcation occurs (SN1 in Fig. 6(b)). When gBK is increased further, the saddle-node 
splits into a folded node (B) and a folded saddle (C) on L + , as shown for gBK = 1 
nS in Fig. 6(c). 

The folded node (B) and the saddle point (A) coalesce at a transcritical bifurcation 
(type II folded saddle-node) when gBK = 3.96 nS (TR1 in Fig. 6(d)). Beyond this, 
the ordinary singularity (A) is a stable node that lies on the top sheet of the criti- 
cal manifold. When gBK =20 nS the folded singularities are either saddles or foci, 
Fig. 6(e). For gBK ~ 32.12 nS the two folded foci on L~ change to folded nodes. 
Finally, when gBK is increased to 32.1224 nS, the fold curves L + and L~ merge. As 
a result, the folded saddles coalesce with the folded nodes at type I folded saddle- 
node bifurcations (SN3 and SN4 in Fig. 6(f)). Beyond this, there is only a stable 
node (A in Fig. 6(g)). The disappearance of the L + and L~ curves correspond to the 
disappearance of the fold in the critical manifold. 

The two-parameter bifurcation diagram in Fig. 7 summarizes the variations of the 
bifurcations in Fig. 5 and Fig. 6 over a range of gK and gBK values. The curves 
TR1 and TR2 correspond to the transcritical bifurcations (type II folded saddle-node 
bifurcations), and SN1-SN4 correspond to the saddle-node bifurcations (type I folded 
saddle-node bifurcations). At gBK = 32.1224 nS the L + and L~ lines coalesce into 
a single line. This contains the SN3 and SN4 bifurcations, up until SN3 and SN4 
coalesce at a codimension-2 bifurcation (for gK = 83.7122 nS). For large gK, the 
L + /L~ line contains no folded singularities (dashed line). 

For gK and gBK values in regions A, D and E there is only a stable node and the 
full system is in a depolarized (A) or hyperpolarized (D or E) steady state. In the left 
portion of region C there is a folded focus which becomes a folded node in the right 
portion of C. This family of folded singularities is on L~ . In region D there is a folded 
node on L~ for negative values of c. Region B consists of the folded nodes on L + , 
and it is the key region for the existence of mixed mode oscillations, since 8 > 0 for 
much of this region (shown below). 



5 Twisted slow manifolds and secondary canards 

The folded nodes discussed above are important since they yield small oscillations 
(for C m > 0) in all trajectories entering the singular funnel. In this section we explain 
the genesis of those oscillations (for more details, see [22, 23, 28, 32]). 
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Fig. 6 V -nullclines (green), c-nullclines (orange) and ordinary and folded singularities for a range of g B % 
values with g K = 7.588 nS. (a) g BK = 0.2 nS, (b) g BK = 0.4 nS, (c) g BK = 1 nS, (d) g BK = 3.96 nS, 
( e ) gBK — 20 nS, (f) g B K = 32.1224 nS, and (g) g B % = 32.2 nS. The color convention for equilibria is 
as in Fig. 5. 
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Fig. 7 Two-parameter 
bifurcation structure for the 
desingularized system. The 
curves TR1 and TR2 represent 
the transcritical bifurcations 
(type II folded saddle-nodes). 
The curves SN1-SN4 represent 
saddle-node bifurcations (type I 
folded saddle-nodes). The 
horizontal line is where the fold 
curves L + and L~ coalesce. 
A codimension-2 bifurcation 
occurs at the intersection of the 
SN curves. 




g K (nS) 



Folded nodes or saddles are characterized by the ratio of their eigenvalues. Let X\ 
and A2 be the eigenvalues of the folded singularity on the fold curve L + such that 
\X\ | < 1 . Define fi as 

A2 

In region A of Fig. 7, which consists of folded saddles, fi < 0. On the TR1 curve 
li = 0 since k\ = 0. Folded nodes occur in region B, so fi > 0. For C m > 0, but small, 
a trajectory approaching a folded node will oscillate, due to twists in the attracting 
and repelling sheets of the slow manifold. The maximum number of oscillations is 
given by [23, 32] 



(i + l 

ah-i 



(5.2) 



which is the greatest integer less than or equal to . At a point in region B and 
close to the TR1 curve in Fig. 7, fi > 0 but small. Hence, S max is large. Similarly on 
SN1 fi = 0, so in region B and close to the SN1 curve fi > 0 and small, so S max is 
large. Between these curves fi increases and S max declines. This is shown in Fig. 8 
for the case gBK = 0-4 pS. The small value of fi over the full range of gx values 
in (Fig. 8a) suggests the system is close to a folded saddle-node bifurcation, either 
type I (SN1) or type II (TR1). 

The attracting sheets of the critical manifold (S a ) and the repelling middle sheet 
(S r ) come together at the fold curves L + and L~. For C m > 0, Fenichel theory [34] 
tells us that the critical manifold perturbs smoothly to invariant attracting (S a ,c m ) and 
repelling (S r ,c m ) manifolds away from L + and L~. However, the critical manifold 
is non-hyperbolic on L + and L~, and perturbs to twisted sheets near these curves 
to preserve uniqueness of solutions [23, 35]. Figure 9 shows how the top attracting 
S+ c (blue) and middle repelling S r ,c m (red) sheets of the slow manifold intersect 
and twist. The numerical method used to compute the slow manifolds was developed 
by Desroches et al. [36, 37]. 

The primary weak canard corresponds to the weak eigendirection of the folded 
node. It is at the intersection of the invariant manifolds S~J~ r and S r c m and serves 

Springer 



Page 14 of 23 



Teka et al. 




Fig. 8 The effects of varying g% on the eigenvalue ratio (fi) and the maximum number of oscillations 
(Smax) for gBK = 0-4 nS. (a) \i = 0 at TR1 and SN1. (b) S max is largest near the bifurcation points. 




Fig. 9 A portion of the twisted slow manifold for C m =2 pF, g% = 4 nS and gBK = 0.4 nS. The top 
attracting (S^~ c , blue surface) and middle repelling (S r> c m > re d surface) sheets of the slow manifold 
are twisted around the blue dashed curve, which is the axis of rotation. The primary strong canard (SC, 
green) moves from the attracting to the repelling sheet without any rotations. The secondary canards £i 
(gray curve, one rotation), £2 (purple curve, two rotations) and £3 (gold curve, three rotations) flow from 
the attracting to repelling sheet with different numbers of rotations. A portion of the pseudo-plateau burst 
trajectory (PPB, black curve) is superimposed and has two small oscillations. The full system has unstable 
equilibrium (cyan, filled curcle). 



as their axis of rotation. All other canards twist about the primary weak canard; they 
follow S+ c as it twists and then follow S r ,c m for a distance as it twists. The primary 
strong canard, which corresponds to the strong eigendirection of the folded node, 
moves along S^" c to S r ,c m without any rotation (SC, green curve in Fig. 9). Other, 
secondary, canards rotate a number of times, depending on how close they are to 
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Fig. 10 (a) Pseudo-plateau burst trajectories are projected onto the (c, V) -plane for g% = 4 nS and 
g BK = 0.4 nS, and different values of C m . Key structures from the desingularized system are also shown. 
WED (pink curve) is the line tangent to the weak eigendirection of the folded node, (b) Magnification of 
panel (a) in the vicinity of the weak eigendirection and fold curve L + . 

the primary strong canard. A secondary canard that makes k small rotations in the 
vicinity of the folded node is called the k th secondary canard. Figure 9 shows the 
first gray), second (§2, purple) and third (§3, olive) secondary canards that make 
one, two and three rotations, respectively. For C m > 0, but small, there are S max — 1 
secondary canards which divide the funnel region between the primary canards into 
^max subsectors [24]. The first subsector is bounded by the strong canard SC and the 
first secondary canard §1 and trajectories entering here have one rotation. The second 
subsector is bounded by §1 and §2 and trajectories entering here have two rotations. 
The last subsector is bounded by the last secondary canard and the primary weak 
canard. The maximal rotation number is achieved in the last subsector; trajectories 
entering here have S max rotations [23, 28, 32]. 

Figure 9 also shows a portion of the pseudo-plateau burst trajectory (PPB, black 
curve) for C m = 2 pF. It enters the funnel region in the rotational subsector bounded 
by §1 and §2, and hence, makes two full rotations and then leaves the repelling sheet 
as it moves towards the lower attracting manifold S~ r . These rotations are the small 
oscillations or "spikes" during the burst active phase. 

Figure 10(a) shows burst trajectories for three values of C m projected onto the 
(c, V) -plane. Also shown are L + , L~, the singular strong canard SC and the folded 
node of the desingularized system. Finally, the line along the weak eigendirection of 
the folded node is included (WED, pink curve). With C m = 0.001 pF the system is 
nearly singular and the "bursting" trajectory enters and leaves the folded node along 
the WED. The small oscillations near the folded node are too small to see. The region 
near the folded node is magnified in Fig. 10(b). With C m = 0. 1 pF the burst trajectory 
again passes through the folded node along the WED, but now the small oscillations 
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g K (nS) g K (nS) 
(a) (b) 

Fig. 11 (a) Mixed mode oscillation borders for C m -> 0. The region of mixed mode oscillations (MMOs) 
is bounded by the two curves TR1 and 8 = 0. Steady state and spiking solutions occur to the left of the 
TR1 curve and to the right of the 8 = 0 curve, respectively, (b) Magnification of panel (a) for a smaller 
range of values of g% and gBK- 



are visible in Fig. 10(b). The small oscillations of this burst trajectory first decrease 
and then increase in amplitude. This is often seen in mixed mode oscillations that 
are associated with a folded node singularity, in contrast to those associated with a 
singular Hopf bifurcation, where the amplitude of successive small oscillations in- 
creases [38]. Finally, with C m = 2 pF (the value used in Fig. 9) the small oscillations 
are prominent even in the larger vertical scale used in Fig. 10(a). 



6 The boundaries of mixed mode oscillations 

For a periodic mixed mode oscillation (i.e., pseudo-plateau bursting) solution to exist, 
there must be a folded node singularity and the periodic orbit must enter the funnel. 
In this section we construct curves in the two-parameter gK~gBK plane that form 
boundaries for the existence of mixed mode oscillations. 

From Fig. 7 we know that folded node singularities only occur in regions B and C 
(and in region D for negative values of the Ca 2+ concentration). Those in region C 
occur on L - and the periodic orbit never enters the corresponding singular funnel. We 
therefore focus on region B. This region is highlighted in Fig. 11(a). Above the TR1 
curve the system has a depolarized stable steady state. Below the SN1 curve the sys- 
tem spikes continuously. Between these curves a folded node singularity exists, and 
the requirement for periodic mixed mode oscillations is that 8 > 0. That is, the sin- 
gular orbit must enter the singular funnel. Thus, the final curve delimiting the MMOs 
region is 8 = 0 (the set of gx and gBK values at which the singular periodic orbit 
intersects both the strong canard and the curve P(L - )), shown in green in Fig. 11. 
For parameter values between the 8 = 0 and TR1 curves periodic mixed mode os- 
cillations, i.e., pseudo-plateau bursting, exist and are stable. This critical region is 
magnified in Fig. 11(b). 

Figure 12 shows how the burst duration and the number of spikes in a burst vary 
over a range of gK and gBK values for C m = 5 pF. A similar map of parameter 
space was used previously in the analysis of a parabolic burster [39]. Two-parameter 
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Fig. 12 The active phase duration and the number of spikes per burst of the full system for C m = 5 pF. 
The system displays steady state, spiking or bursting solutions. Steady state and spiking solutions are 
represented by black dots and small black circles, respectively. The bursting region is bounded by the 
supercritical Hopf bifurcation (HB, black) and the right branch of the period doubling (PD, green) curves. 
The bursting patterns in this region are represented by colored circles. The size of a circle represents the 
active phase duration, with larger circles corresponding to longer active phase durations. The color of a 
circle represents the number of spikes per burst. Cyan circles correspond to smaller number of spikes per 
burst (minimum of two spikes) and the largest dark red circle corresponds to the largest number of spike 
per burst (36 spikes in a burst). 

bifurcation curves of the full system (Eqs. (2.1)-(2.3)) are also shown. These include 
a curve of supercritical Hopf bifurcations (HB, black) and a curve of period doublings 
(PD, green). To the left of the HB curve the system is at a steady state (black dots), and 
to the right of and above the PD curve the system produces continues spiking (small 
black circles). For the values of gK and gBK inside the PD curve the system produces 
pseudo-plateau bursting oscillations (MMOs), represented by colored circles. 

In the bursting region the active phase duration and the number of spike per burst 
vary with respect to the values of gK and gBK ■ The size of each circle represents the 
active phase duration, and the color of the circle (from cyan to dark red) represents 
the number of spikes in a burst. A burst with larger number of spikes has longer 
active phase duration, and in an actual cell this determines the amount of Ca 2+ influx 
and hormone released. The bursts that have the shorter active phase duration and the 
smaller number of spikes occur near the right branch of the PD curve. These bursts 
are represented by smaller cyan circles in Fig. 12. For example, when gBK = 1 nS 
and gK = 6 nS the system produces bursting oscillations with three spikes per burst 
(as in Fig. 4(f)). When one moves away from the right to the left branch of the PD 
curve by increasing gBK ox decreasing gK the burst duration becomes longer and 
the number of spikes in a burst becomes larger. The longest active phase duration is 
about 8.4 sec and the largest number of spikes per burst is about 36, represented by 
the largest dark red circle. These values will change when C m is changed. 

The region between the HB and the left branch of the PD curves is bistable be- 
tween bursting and continuous spiking. Orange circles with small black circles at the 
centers represent bistable solutions that are simulated by varying the initial condi- 
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tions. This shows that the borders of the bursting region are the HB and the right 
branch of the PD curves. The dark blue circles represent bursting oscillations without 
small oscillations since the amplitudes of the spikes are almost zero, i.e., the small 
oscillations are too small to see. 

The results that are shown in Fig. 12 are very consistent with the analysis of the 
mixed mode oscillations in Fig. 11. The HB and TR1 curves overlap, demonstrating 
that for small C m the HB of the full system corresponds to a type II saddle-node 
bifurcation of the desingularized system. Also, the HB curve and the left branch of 
the PD curve are almost indistinguishable for small C m . For these C m values (C m < 
0.001 pF), the right branch of the PD curve converges to the 8 = 0 curve of the 
desingularized system. Hence, the left and right borders of the MMOs in the singular 
limit C m — >► 0 pF correspond to the left and right borders of the bursting region of 
the full system for C m > 0, with the exception that the bursting region is smaller for 
larger values of C m . Also, the bistable region between the PD and HB curves only 
exists as the left PD moves away from the HB, which occurs as C m is increased. 

In Fig. 1 1 the MMOs region delimited by the TR1 and 8 = 0 curves can be divided 
into subregions that have different numbers of small oscillations. For parameter val- 
ues in the subregion near the curve 8 = 0 the periodic orbit enters the funnel region 
near the strong (primary) canard. This subregion corresponds to the first subsector 
of the funnel region, and for C m > 0 only one small oscillation occurs in a burst. 
This corresponds to the jump from the lower attracting sheet to the upper attracting 
sheet and is not due to the folded node. When one moves leftward by decreasing gK , 
8 increases and the periodic orbit enters the funnel region through other subsectors. 
As a result, the number of small oscillations in a burst increases. When one moves 
to the subregion near or on the TR1 curve by decreasing gK further, the periodic 
orbit enters the funnel region through the last subsector. The number of small oscil- 
lations is closer to S max , the maximum number of spikes in a burst as determined 
by the eigenvalues of the folded node. Moreover, increasing gBK has the same effect 
as decreasing gK • These trends in the number of small oscillations obtained from an 
analysis of the desingularized system [28] are expressed far from the singular limit as 
shown in Fig. 12 where C m = 5 pF. Here the longest bursts occur near the HB curves, 
as predicted. 

7 A comparison with a two-fast/one-slow variable analysis 

Using a one-fast/two- slow variable analysis we have shown the genesis of the spikes 
in a burst and how the number of spikes in a burst varies in the gK~gBK parameter 
space. The regions for steady states, pseudo-plateau bursting (mixed mode oscilla- 
tions) and spiking are clearly identified in this parameter space (Fig. 11). This has 
been done by investigating the qualitative changes of the desingularized system when 
parameters gK (Fig. 5) and gBK (Fig- 6) are varied, which are summarized in Fig. 7. 

Here we investigate whether this information can be obtained from a standard two- 
fast/one-slow variable analysis. Figure 13(a) shows a bifurcation diagram of the V-n 
fast subsystem with c treated as a parameter (referred to as a "z-curve"). The subsys- 
tem is bistable over a large range of c values, with stable depolarized and hyperpolar- 
ized steady states, separated by saddle points. The c-nullcline is superimposed, now 
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Fig. 13 Two-fast/one-slow analysis for g^K — 0-4 nS, C m = 10 pF and different values of g% . The black 
"z-curve" is the curve of equilibria of the V-n fast subsystem. This has stable (solid) and unstable (dashed) 
branches, (a) g% =0.1 nS, the full system with stable equilibrium (A\) is in a depolarized steady state, 
(b) gg = 4 nS, the fast subsystem has an unstable limit cycle that emerges from the subcritical Hopf bifur- 
cation (subHB). Pseudo-plateau bursting (PPB, black trajectory) is produced. The equilibrium of the full 
system (A3) is unstable, (c) g% = 5.1 nS, the full system produces periodic spiking that appears unrelated 
to the fast subsystem bifurcation structure. The equilibrium point of the full system (A) is unstable. In all 
panels the system is bistable over a range of c- values. The points A\ and A3 are the same equilibrium 
points as Ai and A3 in Fig. 5, respectively. 



thinking of c as a slowly-changing variable rather than as a parameter. This is the 
standard approach used in a two-fast/one- slow variable analysis. In all three panels 
of Fig. 13 parameters are set at gBK = 0-4 nS, C m = 10 pF, and gK is varied. 

In Fig. 13(a), with gK =0.1 nS, there is an intersection of the c-nullcline on 
the upper stable branch at location Ai. This is a stable equilibrium of the full 3- 
dimensional system, and corresponds to Ai in the analysis shown in Fig. 5. Thus, 
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both types of analysis indicate that the system will come to rest at a depolarized 
steady state when gK =0.1 nS. 

When gK is increased there is a subcritical Hopf bifurcation on the upper branch 
with emergent unstable periodic solutions of the fast- subsystem. This is shown in 
Fig. 13(b) for the case gK = 4 nS. Pseudo-plateau bursting occurs for this and nearby 
values of gK - The full system unstable equilibrium (A3) corresponds to A3 in Fig. 5. 

The superimposed burst trajectory in Fig. 13(b) only weakly follows the fast- 
subsystem bifurcation diagram. Most notably, there are no stable periodic solutions 
of the fast subsystem, only bistability between two steady states. Also, the trajec- 
tory never follows the lower branch of stationary solutions and greatly overshoots the 
lower knee. 

The subcritical Hopf bifurcation migrates leftward when gK is increased to 5.1 nS. 
The unstable branch of periodics goes through a saddle-node bifurcation, yielding a 
branch of stable periodic solutions of the fast subsystem (Fig. 13(c)). There is bista- 
bility between upper and lower branches of the z-curve which is typically a necessary 
condition for bursting with this type of analysis. However, bursting is not produced 
for this value of gK • Instead, the system spikes continuously. 

This example illustrates that features well described by the one-fast/two- slow vari- 
able analysis are not at all well described by a standard two-fast/one- slow variable 
analysis. Most notably, the transition from bursting to spiking is well characterized in 
the one-fast/two-slow variable analysis as the point at which 8 = 0. Note that this is 
not a bifurcation point of the desingularized system, but reflects the jump point from 
the lower sheet of the slow manifold to the upper sheet. In contrast, the bursting to 
spiking transition is not predicted from the two-fast/one-slow analysis, and indeed the 
periodic spiking trajectory of the full system occurs over a range of the fast-subsystem 
bifurcation diagram that contains only stable equilibria. The one-fast/two-slow ap- 
proximation is good even at higher values of C m , for example, when C m = 5 pF 
(Fig. 12). Similar remarks apply for smaller values of C m , where the one-fast/two- 
slow approximation becomes more accurate while the two-fast/one- slow approxima- 
tion does not. The two-fast/one-slow approximation becomes more accurate when c 
is much slower than both V and n, but in this case only a stable steady solution or a 
relaxation oscillation is produced. 



8 Discussion 

The canard mechanism has been used to understand mixed mode oscillations in sev- 
eral neuronal models [30, 37, 40-44]. In these examples, the small oscillations cor- 
respond to subthreshold oscillations that occur between the electrical impulses. We 
have previously analyzed pseudo-plateau bursting in a pituitary lactotroph model us- 
ing canard theory [28]. However, the model used was a simplification in which the 
cytosolic free Ca 2+ concentration was treated as a fixed parameter and the second 
slow variable (in addition to the variable n used here) was an inactivation variable 
for an A-type K + current. In the current paper, we again focused on pseudo-plateau 
bursting in a pituitary lactotroph model, but now with emphasis on a BK-type K + 
current. In this analysis, we have examined the effects of changing the parameters 
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C m , gK and gBK- The parameter gBK is important for producing bursting oscilla- 
tions in actual pituitary cells in which bursting is converted to spiking when BK-type 
K + channels are blocked [45]. 

Here, using C m to control the separation in time scales, we identified two slow 
variables (n, c) and one fast variable (V). Using the one-fast/two- slow variable anal- 
ysis we showed that pseudo-plateau bursting is a canard-induced mixed mode oscilla- 
tion. There are two main requirements for the existence of these bursting oscillations 
[22-24, 32]. One is that the desingularized system must have a folded node singular- 
ity, i.e., the eigenvalue ratio (fi) has to be positive. The second requirement is that the 
singular periodic orbit should enter the singular funnel and pass through the folded 
node, i.e., 8 should be positive. In short, canard-induced mixed mode oscillations 
exist if both /x and 8 are positive. 

Using this technique we can understand several features of the burst and several 
trends that occur as parameters are varied. When both /x and 8 are positive, small 
oscillations are produced during the active phase of a burst and their amplitude is 
proportional to ^/C^ for C m sufficiently small [23]. We obtained the bursting borders 
in the (gK , gBK) -plane (Figs. 1 1 and 12), and predicted how the active phase duration 
and the number of spikes per burst vary with changes in parameters. 

The singular perturbation analysis performed here is technically more effective 
and informative in the singular limit (i.e., for sufficiently small values of C m ) [22, 
23]. However, it provides useful information even far from this limit, as we showed 
in Figs. 11 and 12. Eventually, as the singular parameter (C m ) is increased suffi- 
ciently, new dynamics will be introduced, and the insights from the singular analysis 
are no longer valid. 

The one-fast/two- slow decomposition used here contrasts with the two-fast/one- 
slow variable analysis used previously for pseudo-plateau bursting [10-13]. Our anal- 
ysis explains the origin of the small-amplitude spikes that occur during the active 
phase of pseudo-plateau bursting, the transition between spiking and bursting, and 
information about how the number of spikes per burst varies with parameters. While 
the two-fast/one- slow variable analysis provides little information on these things, 
it does provide valuable information about how one can make a transition between 
plateau and pseudo-plateau bursting as one or more parameters are changed [12]. It 
also provides information about complex phase resetting properties [11] and the ter- 
mination of spikes in a burst [17]. Both fast/slow decompositions are approximations, 
however, to a system that evolves on three time scales. Some studies [13, 17, 18] fo- 
cus on the dynamics of the full system, and illustrate the complexity of the seemingly 
simple set of equations. The advantage of obtaining useful information of the full sys- 
tem by a two-fast/one- slow or one-fast/two-slow decomposition points to the fact that 
system (2.1)-(2.3) actually evolves on three time scales: V fast, n intermediate and c 
slow. This can also be seen by the magnitude of fi which is bounded from above by 
Umax ~ 0.07 (Fig. 8(a)). Hence, we are close to folded saddle-node regimes (type I 
and type II) [33,38] and a more detailed bifurcation analysis may explain the relation 
between the two-fast/one-slow and one-fast/two- slow splitting. This is left for future 
work. 



Springer 



Page 22 of 23 



Teka et al. 



Competing interests 

The authors declare that they have no competing interests. 



Authors' contributions 

WT, JT, TV, MW, and RB performed the analysis. WT wrote the manuscript with assistance from JT, TV, 
MW, and RB. All authors read and approved the final manuscript. 

Acknowledgements This work was supported by NSF grant DMS 0917664 to RB and NIH grant DK 
043200 to RB and JT. 



References 

1. Lisman, J.E.: Bursts as a unit of neural information: making unreliable synapses reliable. Trends 
Neurosci. 20, 38-43 (1997) 

2. Van Goor, R, Zivadinovic, D., Martinez-Fuentes, A.J., Stojilkovic, S.S.: Dependence of pituitary hor- 
mone secretion on the pattern of spontaneous voltage-gated calcium influx. Cell type-specific action 
potential secretion coupling. J. Biol. Chem. 276, 33840-33846 (2001) 

3. Stojilkovic, S.S., Zemkova, H., Van Goor, F.: Biophysical basis of pituitary cell type-specific Ca 2+ 
signaling-secretion coupling. Trends Endocrinol. Metab. 16, 152-159 (2005) 

4. Kuryshev, Y.A., Childs, G.V, Ritchie, A.K.: Corticotropin-releasing hormone stimulates Ca 2+ entry 
through L- and P-type Ca 2+ channels in rat corticotropes. Endocrinology 137, 2269-2277 (1996) 

5. Tsaneva-Atanasova, K., Sherman, A., Van Goor, E, Stojilkovic, S.S.: Mechanism of spontaneous and 
receptor-controlled electrical activity in pituitary somatotrophs: experiments and theory. J. Neuro- 
physiol. 98, 131-144(2007) 

6. Falke, L.C., Gillis, K.D., Pressel, D.M., Misler, S.: Perforated patch recording allows long-term mon- 
itoring of metabolite-induced electrical activity and voltage-dependent Ca 2+ currents in pancreatic 
islet p cells. FEBS Lett. 251, 167-172 (1989) 

7. Kinard, T., Vries, G.D., Sherman, A., Satin, L.S.: Modulation of the bursting properties of single 
mouse pancreatic beta-cells by artificial conductances. Biophys. J. 76(3), 1423-1435 (1999) 

8. Zhang, M., Goforth, P., Bertram, R., Sherman, A., Satin, L.: The Ca 2+ dynamics of isolated mouse 
/3-cells and islets: implications for mathematical models. Biophys. J. 84, 2852-2870 (2003) 

9. LeBeau, A.P., Robson, A.B., McKinnon, A.E., Sneyd, J.: Analysis of a reduced model of corticotroph 
action potentials. J. Theor. Biol. 192, 319-339 (1998) 

10. Tabak, J., Toporikova, N., Freeman, M.E., Bertram, R.: Low dose of dopamine may stimulate prolactin 
secretion by increasing fast potassium currents. J. Comput. Neurosci. 22, 211-222 (2007) 

11. Stern, J.V, Osinga, H.M., LeBeau, A., Sherman, A.: Resetting behavior in a model of bursting in 
secretory pituitary cells: distinguishing plateaus from pseudo-plateaus. Bull. Math. Biol. 70, 68-88 
(2008) 

12. Teka, W., Tsaneva-Atanasova, K., Bertram, R., Tabak, J.: From plateau to pseudo-plateau bursting: 
making the transition. Bull. Math. Biol. 73(6), 1292-1311 (2011) 

13. Tsaneva-Atanasova, K., Osinga, H.M., Riess, T., Sherman, A.: Full system bifurcation analysis of 
endocrine bursting models. J. Theor. Biol. 264, 1133-1146 (2010) 

14. Safiulina, V.F., Zacchi, P., Taglialatela, M., Yaari, Y., Cherubini, E.: Low expression of Kv7/M chan- 
nels facilitates intrinsic and network bursting in the developing rat hippocampus. J. Physiol. 586(22), 
5437-5453 (2008) 

15. Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: Teramoto, E., 
Yamaguti, M. (eds.) Mathematical Topics in Population Biology, Morphogenesis, and Neurosciences. 
Lecture Notes in Biomathematics, pp. 267-281. Springer, Berlin (1987) 

16. Rinzel, J., Ermentrout, G.B.: Analysis of neural excitability and oscillations. In: Koch, C., Segev, I. 
(eds.) Methods in Neuronal Modeling: From Synapses to Networks, 2nd edn, pp. 251-292. MIT Press, 
Cambridge, MA (1998) 



4i) Springer 



Journal of Mathematical Neuroscience (2011) 1:12 



Page 23 of 23 



17. Nowacki, J., Mazlan, S., Osinga, H.M., Tsaneva-Atanasova, K.: The role of large-conductance 
calcium-activated K + (BK) channels in shaping bursting oscillations of a somatotroph cell model. 
Physica D 239, 485-493 (2010) 

18. Osinga, H.M., Tsaneva-Atanasova, K.: Dynamics of plateau bursting depending on the location of its 
equilibrium. J. Neuroendocrinol. 22(12), 1301-1314 (2010) 

19. Bertram, R., Sherman, A.: Negative calcium feedback: the road from Chay-Keizer. In: Coombes, S., 
Bressloff, P. (eds.) The Genesis of Rhythm in the Nervous System, pp. 19-48. World Scientific, New 
Jersey (2005) 

20. Butera, R.J., Rinzel, J., Smith, J.C.: Models of respiratory rhythm generation in the pre-Botzinger 
complex. I. Bursting pacemaker neurons. J. Neurophysiol. 82, 382-397 (1999) 

21. Del Negro, C.A., Hsiao, C.F., Chandler, S.H.: Outward currents influencing bursting dynamics in 
guinea pig trigeminal motoneurons. J. Neurophysiol. 81(4), 1478-1485 (1999) 

22. Szmolyan, P., Wechselberger, M.: Canards in R 3 . J. Differ. Equ. 177, 419-453 (2001) 

23. Wechselberger, M.: Existence and bifurcation of canards in R 3 in the case of a folded node. SIAM J. 
Appl. Dyn. Syst. 4, 101-139 (2005) 

24. Br0ns, M., Krupa, M., Wechselberger, M.: Mixed mode oscillations due to the generalized canard 
phenomenon. Fields Inst. Commun. 49, 39-63 (2006) 

25. Desroches, M., Guckenheimer, J., Krauskopf, B., Kuehn, C, Osinga, H., Wechselberger, M.: Mixed- 
mode oscillatons with multiple time-scales. SIAM Rev. (201 1, in press) 

26. Ermentrout, B.: Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for 
Researchers and Students. SIAM, Philadelphia (2002) 

27. Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Oldeman, K.E., Paffenroth, R.C., 
Sandstede, B., Wang, X.J., Zhang, C: AUTO-07P:: Continuation and bifurcation software for ordinary 
differential equations, http://cmvl.cs.concordia.ca/ (2007) 

28. Vo, T., Bertram, R., Tabak, J., Wechselberger, M.: Mixed mode oscillations as a mechanism for 
pseudo-plateau bursting. J. Comput. Neurosci. 28, 443-458 (2010) 

29. Wechselberger, M.: A propos de canards (apropos canards). American Mathematical Society (2011, 
in press) 

30. Rubin, J., Wechselberger, M.: Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley 
model. Biol. Cybern. 97, 5-32 (2007) 

31. Szmolyan, P., Wechselberger, M.: Relaxation oscillations in R 3 . J. Differ. Equ. 200, 69-104 (2004) 

32. Rubin, J., Wechselberger, M.: The selection of mixed-mode oscillations in a Hodgkin-Huxley model 
with multiple timescales. Chaos 18, 015105 (2008) 

33. Krupa, M., Wechselberger, M.: Local analysis near a folded saddle-node singularity. J. Differ. Equ. 
248(12), 2841-2888 (2010) 

34. Fenichel, N.: Geometric singular perturbation theory. J. Differ. Equ. 31, 53-98 (1979) 

35. Guckenheimer, J., Haiduc, R.: Canards at folded nodes. Mosc. Math. J. 5, 91-103 (2005) 

36. Desroches, M., Krauskopf, B., Osinga, H.M.: The geometry of slow manifolds near a folded node. 
SIAM J. Appl. Dyn. Syst. 7(4), 1131-1162 (2008) 

37. Desroches, M., Krauskopf, B., Osinga, H.M.: Mixed-mode oscillations and slow manifolds in the 
self-coupled FitzHugh-Nagumo system. Chaos 18, 015107 (2008) 

38. Guckenheimer, J.: Singular Hopf bifurcation in systems with two slow variables. SIAM J. Appl. Dyn. 
Syst. 7(4), 1355-1377 (2008) 

39. Guckenheimer, J., Gueron, S., Harriswarrick, R.M.: Mapping the dynamics of a bursting neuron. 
Philos. Trans. R. Soc. Lond. B, Biol. Sci. 341(1298), 345-359 (1993) 

40. Guckenheimer, J., HarrisWarrick, R., Peck, J., Willms, A.: Bifurcation, bursting, and spike frequency 
adaptation. J. Comput. Neurosci. 4(3), 257-277 (1997) 

41. Br0ns, M., Kaper, T.J., Rotstein, H.G.: Introduction to focus issue: mixed mode oscillations: experi- 
ment, computation, and analysis. Chaos 18, 015101 (2008) 

42. Erchova, I., McGonigle, D.J.: Rhythms of the brain: An examination of mixed mode oscillation ap- 
proaches to the analysis of neurophysiological data. Chaos 18, 015115 (2008) 

43. Rotstein, H.G., Wechselberger, M., Kopell, N.: Canard Induced mixed-mode oscillations in a medial 
entorhinal cortex Layer II Stellate Cell Model. SIAM J. Appl. Dyn. Syst. 7(4), 1582-1611 (2008) 

44. Ermentrout, B., Wechselberger, M.: Canards, clusters, and synchronization in a weakly coupled in- 
terneuron model. SIAM J. Appl. Dyn. Syst. 8, 253-278 (2009) 

45. Van Goor, F, Li, Y.X., Stojilkovic, S.S.: Paradoxical role of large-conductance calcium-activated K + 
(BK) channels in controlling action potential-driven Ca 2+ entry in anterior pituitary cells. J. Neurosci. 
21, 5902-5915 (2001) 



Springer 



